#include <cmath>

#include "../Iteration_Params.hxx"

Okada::Displacement Iteration_Params::ua
(const FTensor::Tensor1<double,3> &dislocation,
 const double &alp1, const double &alp2,
 const double &sin_dip, const double &cos_dip) const
{
  Okada::Displacement result;
  FTensor::Index<'a',3> a;
  FTensor::Index<'b',3> b;
  result.first(a)=0;
  result.second(a,b)=0;

  double xy(xi*y11), qx(q*x11), qy(q*y11);

  const double pi(4*atan(1.0));
  FTensor::Tensor1<double,3> first;
  FTensor::Tensor2<double,3,3> second;
  /* STRIKE-SLIP CONTRIBUTION */
  if(dislocation(0)!=0)
    {
      first(0)=tt/2 + alp2*xi*qy;
      first(1)=       alp2*q/r;
      first(2)=alp1*ale - alp2*q*qy;

      second(0,0)=-alp1*qy  -alp2*xi2*q*y32;
      second(0,1)=          -alp2*xi*q/r3;
      second(0,2)= alp1*xy  +alp2*xi*q2*y32;
      second(1,0)= alp1*xy*sin_dip        +alp2*xi*fy+d/2*x11;
      second(1,1)=                    alp2*ey;
      second(1,2)= alp1*(cos_dip/r+qy*sin_dip) -alp2*q*fy;
      second(2,0)= alp1*xy*cos_dip        +alp2*xi*fz+y/2*x11;
      second(2,1)=                    alp2*ez;
      second(2,2)=-alp1*(sin_dip/r-qy*cos_dip) -alp2*q*fz;

      result.first(a)=first(a)*dislocation(0)/(2*pi);
      result.second(a,b)=second(a,b)*dislocation(0)/(2*pi);
    }

  /* DIP-SLIP CONTRIBUTION */
  if(dislocation(1)!=0)
    {
      first(0)=           alp2*q/r;
      first(1)=    tt/2 +alp2*et*qx;
      first(2)= alp1*alx -alp2*q*qx;
      second(0,0)=        -alp2*xi*q/r3;
      second(0,1)= -qy/2 -alp2*et*q/r3;
      second(0,2)= alp1/r +alp2*q2/r3;
      second(1,0)=                      alp2*ey;
      second(1,1)= alp1*d*x11+xy/2*sin_dip +alp2*et*gy;
      second(1,2)= alp1*y*x11          -alp2*q*gy;
      second(2,0)=                      alp2*ez;
      second(2,1)= alp1*y*x11+xy/2*cos_dip +alp2*et*gz;
      second(2,2)=-alp1*d*x11          -alp2*q*gz;

      result.first(a)+=first(a)*dislocation(1)/(2*pi);
      result.second(a,b)+=second(a,b)*dislocation(1)/(2*pi);
    }

  /* TENSILE-FAULT CONTRIBUTION */
  if(dislocation(2)!=0)
    {
      first(0)=-alp1*ale -alp2*q*qy;
      first(1)=-alp1*alx -alp2*q*qx;
      first(2)=    tt/2 -alp2*(et*qx+xi*qy);
      second(0,0)=-alp1*xy  +alp2*xi*q2*y32;
      second(0,1)=-alp1/r   +alp2*q2/r3;
      second(0,2)=-alp1*qy  -alp2*q*q2*y32;
      second(1,0)=-alp1*(cos_dip/r+qy*sin_dip)  -alp2*q*fy;
      second(1,1)=-alp1*y*x11         -alp2*q*gy;
      second(1,2)= alp1*(d*x11+xy*sin_dip) +alp2*q*hy;
      second(2,0)= alp1*(sin_dip/r-qy*cos_dip)  -alp2*q*fz;
      second(2,1)= alp1*d*x11         -alp2*q*gz;
      second(2,2)= alp1*(y*x11+xy*cos_dip) +alp2*q*hz;

      result.first(a)+=first(a)*dislocation(2)/(2*pi);
      result.second(a,b)+=second(a,b)*dislocation(2)/(2*pi);
    }
  return result;
}
